TK!2
=v
#1
:n
L1
:v
10
:s
i
:c
spring length
#2
:n
L2
:v
10
:s
i
:c
   "     "
#3
:n
K1
:v
8
:s
i
:c
spring stiffness
#4
:n
K2
:v
1
:s
i
:c
   "      "
#5
:n
P1
:v
5
:s
i
:c
load in x1 direction
#6
:n
P2
:v
5
:s
i
:c
load in x2 direction
#8
:n
n
:v
2
:s
i
:c
number of unknowns
#9
:n
OFN
:v
'springs
:s
i
:c
objective function name
#10
:n
dOFN
:v
'dsprings
:s
i
:c
function returning derivatives of OFN
#11
:n
x10
:v
3
:s
i
:c
initial guess for x1
#12
:n
x20
:v
3
:s
i
:c
initial guess for x2
#13
:n
TOL
:v
.000001
:s
i
:c
convergence tolerance
=u
=l
#1
:n
x
:v
8.632508083036267
4.53256322250374
E
#2
:n
dx
:v
-.0004092997232647991
.0003797550199304908
E
#3
:n
g
:v
-.00447315910232593
.001889633021635362
E
#4
:n
h
:v
-.0037870482496099
.00351368569721196
E
#5
:n
x1
:v
E
#6
:n
x2
:v
E
#7
:n
PE
:v
E
#8
:n
k
:v
164
E
#9
:n
Item
:v
'Function_calls
'Direction_changes

'Coordinate_x1
'Coordinate_x1

'Object_function_PE
E
#10
:n
Value
:v
E
=f
#1
:n
springs
:c
Potential energy of a two-spring system
:t
Rule
:p
L1,L2,K1,K2,P1,P2
:a
x1,x2
:v
PE
#2
:n
dsprings
:c
Derivatives dPE/dx1 and dPE/dx2
:t
Rule
:p
L1,L2,K1,K2,P1,P2
:a
x1,x2
:v
PE
#3
:n
congrad
:c
Multidim. optimization: conjugate gradient method
:t
Procedure
:p
n,dOFN
:i
TOL
#4
:n
linemin
:c
Line minimization
:t
Procedure
:p
n
:o
fxmin
#5
:n
object
:c
Objective function evaluator
:t
Procedure
:p
OFN
:i
x
:o
y
#6
:n
bracket
:c
Bracketing a minimum
:t
Procedure
:i
FUN,X1,X2
:o
x1,x2,x3,f1,f2,f3
#7
:n
Brent
:c
Inverse parabolic interpolation (Brent's method)
:t
Procedure
:i
FUN,a,b,c,tol
:o
x,fx
#8
:n
trace
:c
Optimization trace
:t
Procedure
:i
x1,x2,PE,Flag
#1
:n
springs
:z
pe1 = K1*(sqrt(x1^2+(L1-x2)^2) - L1)^2
pe2 = K2*(sqrt(x1^2+(L2+x2)^2) - L2)^2
PE = (pe1+pe2)/2 - P1*x1 - P2*x2

" Note: This function returns the potential energy of a two-spring system.
"       The springs with lengths L1 and L2 and stiffnesses K1 and K2 are
"       loaded with forces P1 and P2 along the x- and y-axes; the potential
"       energy is a function of the displacements x1 and x2.

"                                 P2
"                                ³
"                                ÀÄÄÄÄ P1
"                  ±oÄ@@@@@@@@@@ÄoÄ@@@@@@@@@@Äo±

call trace(x1,x2,PE,'record)       " optimization monitor
place('Value,1) = 'Value[1] + 1    " optimization trace
place('Value,4) = x1
place('Value,5) = x2
place('Value,7) = PE
/E
E
#2
:n
dsprings
:z
" Numeric differentiation:
dx = 1e-5
PE = springs(x1,x2)
place('dx,1) = (springs(x1*(1+dx),x2)-PE)/(x1*dx)
place('dx,2) = (springs(x1,x2*(1+dx))-PE)/(x2*dx)

" Symbolic differentiation:
aux1 = sqrt(x1^2+(L1-x2)^2) 
pe1 = K1*(aux1 - L1)^2
aux2 = sqrt(x1^2+(L2+x2)^2) 
pe2 = K2*(aux2 - L2)^2
PE = (pe1+pe2)/2 - P1*x1 - P2*x2

dpe1%daux1 = K1*2*(aux1 - L1)
daux1%dx1 = x1/aux1
daux1%dx2 = -(L1-x2)/aux1
dpe1%dx1 = dpe1%daux1 * daux1%dx1
dpe1%dx2 = dpe1%daux1 * daux1%dx2

dpe2%daux2 = K2*2*(aux2 - L2)
daux2%dx1 = x1/aux2
daux2%dx2 = (L1+x2)/aux2
dpe2%dx1 = dpe2%daux2 * daux2%dx1
dpe2%dx2 = dpe2%daux2 * daux2%dx2

dPE%dx1 = (dpe1%dx1+dpe2%dx1)/2 - P1
dPE%dx2 = (dpe1%dx2+dpe2%dx2)/2 - P2

L^2 = dPE%dx1^2 + dPE%dx2^2
place('dx,1) = dPE%dx1
place('dx,2) = dPE%dx2
/E
1
2
3
4
E
#3
:n
congrad
:z
" Notation:  n     number of unknowns
"            dOFN  name of function returning the objective
"                  function derivatives
"            TOL   convergence tolerance

" Supporting functions: linemin,bracket,Brent,object

" Description: This function represents the essence of the Polak-Ribiere
" conjugate gradient algorithm.  It computes the gradient, calls the function
" linemin  which performs the optimization along the gradient, and checks
" the conditions of termination of the optimization process.

'Value[1]:= 0      " monitor initialization
'Value[2]:= 0           "          "
itmax:= 200
eps:= 1e-10
V:= apply(dOFN,'x[1],'x[2])
for j=1 to n
   'g[j]:= -'dx[j]
   'h[j]:= 'g[j]
   'dx[j]:= 'h[j]
next j
for its=1 to itmax
   minV:= linemin()
   if 2*abs(minV - V) <= TOL*(abs(V) + abs(V) + eps) return
   V:= apply(dOFN,'x[1],'x[2])
   GG:= 0
   DGG:= 0
   for j=1 to n
      GG:= GG + 'g[j]^2
      DGG:= DGG + 'dx[j]^2
   next j
   if GG=0 return
   GAM:= DGG/GG
   for j=1 to n
      'g[j]:= -'dx[j]
      'h[j]:= 'g[j] + GAM*'h[j]
      'dx[j]:= 'h[j]
   next j
next its
/E
E
#4
:n
linemin
:z
" This function launches local minimization along the conjugate gradient, 
" calling first the function  bracket  and then  Brent .

'Value[2]:= 'Value[2] + 1
call bracket('object,0,.1;x1,x2,x3,f1,f2,f3)
(xmin,fxmin):= Brent('object,x1,x2,x3,1e-6)
for i=1 to n
   'dx[i]:= xmin*'dx[i]
   'x[i]:= 'x[i] + 'dx[i]
next i
/E
E
#5
:n
object
:z
" This function transforms a coordinate  x  along the conjugate gradient
" into the objective function (OFN) parameters  x1, x2, ... , and evaluates
" the objective function.

x1:= 'x[1] + x*'dx[1]
x2:= 'x[2] + x*'dx[2]
y:= apply(OFN,x1,x2)
/E
E
#6
:n
bracket
:z
" Notation:  FUN    name of the objective function
"            X1,X2  initial search interval

" Description: The procedure searches for (and returns) a set of three values
" of abscissa  x ,  x1 < x2 < x3  or  x1 > x2 > x3 , such that the corres- 
" ponding function values  f1 > f2  and  f3 > f2 .  This means that the
" minimum should be somewhere between  x1  and  x3 ,  and that an optimiza-
" tion procedure, such as golden section search 'Golden' or inverse parabolic
" interpolation 'Brent' can converge to it.

x1:= X1
x2:= X2
gold:= 1.618034                    " interval magnification ratio
glimit:= 100                       " maximum allowed magnification
f1:= apply(FUN,x1)
f2:= apply(FUN,x2)
if f1 > f2 goto A
  (x1,x2):= (x2,x1)                " securing downhill move from x1 to x2
  (f1,f2):= (f2,f1)
A: x3:= x2 + gold*(x2-x1)
   f3:= apply(FUN,x3)
B: if f2 < f3 return
   xf1:= (x2-x1)*(f2-f3)           " finding x4 by parabolic extrapolation
   xf3:= (x2-x3)*(f2-f1)
   if f1=f3 then f1:= .999*f1               " preventing division by 0
   x4:= x2 - ((x2-x3)*f3 - (x2-x1)*f1) / (f3-f1)
   xlimit:= x2 + glimit*(x3-x2)             " current step limit
   if (x2-x4)*(x4-x3)<=0 goto E
     f4:= apply(FUN,x4)
     if x4 >= x3 goto C
       (x1,f1):= (x2,f2)                    " minimum between x2 and x3
       (x2,f2):= (x4,f4)
       goto B
     C: if f4 <= f2 goto D
          (x3,f3):= (x4,f4)                 " minimum between x1 and x4
          goto B
        D: x4:= x3 + gold*(x3-x2)
           f4:= apply(FUN,x4)
           goto H
   E: if (x3-x4)*(x4-xlimit) <= 0 goto F
        f4:= apply(FUN,x4)                  " x4 is between x3 and xlimit
        if f4 >= f3 goto H
          (x2,f2):= (x3,f3)
          (x3,f3):= (x4,f4)
          x4:= x3 + gold*(x3-x2)
          f4:= apply(FUN,x4)
          goto H
      F: if (x4-xlimit)*(xlimit-x3) < 0 goto G
           x4:= xlimit
           f4:= apply(FUN,x4)
           goto H
         G: x4:= x3 + gold*(x3-x2)
            f4:= apply(FUN,x4)
H: (x1,f1):= (x2,f2)                        " eliminating x1
   (x2,f2):= (x3,f3)
   (x3,f3):= (x4,f4)
   goto B
/E
E
#7
:n
Brent
:z
" Notation:  FUN    name of the objective function
"            a,b,c  values of  x  bracketing the minimum
"            tol    desired relative precision of  x  giving minimum
"                   FUN (shouldn't be less than 1e-8)

" Description:   s e e   m o d e l   OPTIM1D

" Note: Maximum number of iteration steps is set in the statement labeled E.

gold:= .381966
eps:= 1e-10
(a,b):= (min(a,c),max(a,c))
(v,fv):= (b,apply(FUN,b))
(w,fw):= (v,fv)
(x,fx):= (v,fv)
e:= 0
iter:= 0
LOOP:
iter:= iter + 1
m:= (a+b)/2
tol1:= tol*abs(x) + eps
tol2:= 2*tol1
if abs(x-m) <= (tol2-(b-a)/2) return
if abs(e) <= tol1 goto A

  r:= (x-w)*(fx-fv)                      " trying parabolic fit
  q:= (x-v)*(fx-fw)
  p:= (x-v)*q - (x-w)*r
  q:= 2*(q-r)
  if q > 0 then p:= -p
  q:= abs(q)
  etemp:= e
  e:= d
  if or(abs(p)>=abs(.5*q*etemp),p<=q*(a-x),p>=q*(b-x)) goto A

    d:= p/q                               " taking parabolic step
    u:= x + d
    if or((u-a)<tol2,(b-u)<tol2) then d:= abs(tol1)*sgn(m-x)
    goto B
A: if x < m then  e:= b-x  else  e:= a-x

   d:= gold*e                             " taking golden section step
B: if abs(d) < tol1 then  u:= x + abs(tol1)*sgn(d)  else  u:= x + d

   fu:= apply(FUN,u)                      " function evaluation at new point

   if fu > fx goto C                      " incorporating new point
     if u < x then  b:= x  else  a:= x    " into the search scheme
     (v,fv):= (w,fw)
     (w,fw):= (x,fx)
     (x,fx):= (u,fu)
     goto E
   C: if u < x then  a:= u  else  b:= u
      if and(fu>fw,w<>x) goto D
        (v,fv):= (w,fw)
        (w,fw):= (u,fu)
        goto E
      D: if or(fu<=fv,v=x,v=w) then (v,fv):= (u,fu)
E: if iter<50 goto LOOP
error:= convergence_failure
/E
E
#8
:n
trace
:z
" Notation:  x1     unknown variable value
"            x2     unknown variable value
"            PE     objective function value
"            flag   initialization/recording flag

" Description: This is an example of a procedure for tracing the optimiza-
" tion process.  Here are the steps for incorporating it into your model:

" 1. Select the variables you want to record and make them Input Variables
"    in the  trace  function.
" 2. Add a variable  Flag  as an input variable.
" 3. Select the list names under which you want to record the trace infor-
"    mation.  Blank these lists in the initialization part of this procedure
"    and use them on the left hand side of the assignment statement in the
"    recording part. (We use the same names for the variables and lists
"    in this example.)
" 4. Set up some list element as a record counter.  Set it to zero in the
"    initialization part and increment it in the recording part. (This 
"    example uses k[1].)
" 5. Include an initialization call of the function  trace  before the
"    optimization procedure call on the Rule Sheet of Function Subsheet.
"    All the arguments may be set to zero, e.g.  call trace(0,0,0) .
" 6. Include the recording call(s) in the Function Subsheet evaluating the
"    objective function.  The calls must have as a last argument a symbolic
"    constant  'record .
" 7. set up an Interactive Table for displaying the trace lists.

" Note: You may want to display more information than just the values of 
" the unknowns and the objective function.  Any such information can be sent
" to the  trace  via input variable arguments. 

if Flag='record goto RECORD
'k[1]:= 0                         " initialization
  call blank('x1)
  call blank('x2)
  call blank('PE)
  return
RECORD:
'k[1]:= 'k[1] + 1                 " recording the objective function call
'x1['k[1]]:= x1
'x2['k[1]]:= x2
'PE['k[1]]:= PE
/E
E
=r
#1
:r
" ***  Multidimensional Optimization, Part 4 - Conjugate Gradient Method
:s
C
#3
:r
call trace(0,0,0,0)     " trace initialization
#4
:r
place('x,1) = x10       " initial values
#5
:r
place('x,2) = x20       " of unknowns
#6
:r
call congrad(TOL)
#8
:r
" N o t e s :
:s
C
#10
:r
" This model demonstrates the conjugate gradient method of multidimensional
:s
C
#11
:r
" optimization (Polak-Ribiere version) along the lines suggested in 
:s
C
#12
:r
" "Numerical Recipes" by W.H.Press et al. (Cambridge University Press, 1986).
:s
C
#13
:r
" The problem example (minimization of the potential energy of a two-spring
:s
C
#14
:r
" system; see the function  springs ) is described in the model OPTIM2A. The
:s
C
#15
:r
" same problem is also solved by different methods in the models OPTIM2B and
:s
C
#16
:r
" OPTIM2C.  The contour map of the objective function is included in the
:s
C
#17
:r
" model OPTIM2C.
:s
C
#19
:r
" You should invoke the solution (F9) when in the interactive table  monitor,
:s
C
#20
:r
" where you can follow the progress of the optimization process.  Also
:s
C
#21
:r
" included in this model is the function  trace  which records the iteration
:s
C
#22
:r
" path for later perusal (it may be switched off by cancelling the trace
:s
C
#23
:r
" calls in the Rule Sheet and in the Function Subsheet  springs ).
:s
C
#25
:r
" Gradient methods require evaluation of partial derivatives of the objective
:s
C
#26
:r
" function with respect to the unknown variables.  This is performed in this
:s
C
#27
:r
" model by the function  dsprings .  The model in its current state uses
:s
C
#28
:r
" symbolic derivatives, but the numeric differentiation alternative is also
:s
C
#29
:r
" included in the function (and cancelled).  Numerical differentiation needs
:s
C
#30
:r
" to be used when the symbolic derivatives are hard or impossible to obtain.
:s
C
=p
=t
#1
:n
trace
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Iteration trace"
:u
#1
:l
x1
:w
14
:h
""
#2
:l
x2
:w
14
:h
""
#3
:l
PE
:w
14
:h
""
E
#2
:n
monitor
:s
Screen
:v
Vertical
:f
1
:r
""
:c
""
:t
"Optimization Monitor"
:u
#1
:l
Item
:w
21
:h
""
#2
:l
Value
:w
18
:h
""
E
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
No
:a
Yes
:r
Yes
:h
No
:t
1
:n
No
:d
LPT1
%Tr,f,0,1,0,1,0
